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Abstract 

We  estimate  the  uncertainty  associated  with  the  application  of  the 
stong-shock  methodologj'  for  nuclear  yield  determination  for  the  ide¬ 
alized  case  of  a  fully-coupled  event  and  a  spherically  symmetric  wave 
field.  Time-of-arrival  data  from  a  number  of  U.  S.  nuclear  tests  are  ana¬ 
lyzed  in  order  to  estimate  measurement  uncertainties  and  uncertainties 
in  the  algorithmic  parameters.  These  findings  are  then  employed  to 
estimate  resulting  confidence  limits  for  yield  determinations  produced 
via  this  algorithm. 
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1  Introduction 

1.1  Background 

The  hydrodynamic  yield  algorithm  is  often  used  to  estimate  yields  for  under¬ 
ground  nuclear  events.  For  this  algorithm,  it  is  assumed  that,  under  certain 
conditions,  the  yield  of  an  underground  nuclear  explosion  can  be  reasonably 
well  determined  from  measurements  of  the  propagation  of  the  explosion- 
produced  shock  wave  through  the  ambient  geological  medium.  If  the  explo¬ 
sion  is  spherically  symmetric  and  the  measurements  are  taken  in  a  pressure 
regime  which  far  exceeds  the  yield  strength  of  the  medium  (the  so-called 
hydrodynamic  regime),  the  shock  radius  appears  to  grow  as  a  power-law 
function  of  time.  In  particular,  for  some  range  of  t,  the  position  of  the  shock 
front  is  assumed  to  be  given  by  the  generalized  Taylor  blast-wave  solution:’ 


where  i  is  time  after  explosion,  /?,  is  the  distance  from  the  center  of  the  ex¬ 
plosion,  and  W  is  the  yield  of  the  device.  The  parameters  a  and  b  have  been 
empirically  determined  and  are  assumed  to  be  independent  of  the  material 
medium  and  device  yield  over  a  wide  range  of  geologic  media  and  yields. 

’The  original  derivation  by  Taylor,^’^^  for  shock  waves  propagating  in  air.  gave  6=5. 
since  he  assumed  a  gamma-law  gas  with  constant  7.  For  shocks  in  non-ideal  materials 
such  as  rock  a  more  general  form  is  required. 
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This  equation  has  formed  the  basis  of  the  h>drodynamic-yield  determina¬ 
tion  techniques  used  on  L'.S.  nuclear  events  sir«ce  the  early  1960‘s. 

We  have  analyzed  shock  time-of-arrival  (TOA)  data  from  43  U.  S.  nuclear 
events  with  yields  greater  than  40  kt,  detonated  in  a  number  of  different 
media.  Some  of  the  questions  which  we  have  attempted  to  address  in  this 
analysis  are:  1)  Can  we  reliably  estimate  the  yield  of  a  nuclear  event  from 
TO.A.  data  alone,  using  the  "standard’’  hydrodynamic  yield  algorithm  based 
on  Eqn.  (1)?  2)  If  so,  what  are  the  confidence  limits  which  can  be  placed  on 
such  yield  estimates  and  how  do  these  limits  scale  with  variations  in  yield? 
3)  Are  there  any  significant  variations  in  these  results  for  different  working 
point  media  (e.g.,  basalt,  tuff,  granite)? 

In  this  analysis  we  have,  for  simplicity,  assumed  that  the  geometry  is 
spherically  symmetric  and  that  the  events  are  fully  coupled.  For  relatively 
large  yields  in  reasonably  small  zero  room  volumes,  these  are  evidently  good 
assumptions.  As  the  e.\p!osive  yield  becomes  small  or  as  the  cavity  size 
becomes  large,  nonspherical  effects  or  decoupling  may  become  important. 

1.2  The  strong-shock  algorithm 

Here  we  present  a  brief  description  of  the  hydrodynamic  or  strong-shock 
yield  algorithm  as  it  is  presently  applied.!®!  Figure  1  shows  an  example  of 
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a  typical  {R,J)  data  set  from  a  U.  S.  nuclear  event,  obtained  in  this  case 
using  the  SLIFER^^l  measurement  technique.  The  data  are  presented  in 


Shock  Radius  vs  Time 


Figure  I:  Scaled  shock  radius  vs  scaled  time  data  for  a  typical  nuclear  event. 

scaled  form  where  the  scaling  factor  for  both  R,  and  i  is  the  cube  root  of 
the  official  yield  for  the  event,  but  the  details  of  the  algorithm  are  largely 
unaffected  by  this  scaling.  The  most  difficult  part  of  the  algorithm  is  making 
a  proper  determination  of  the  region  of  data  to  analyze.  Equation  (1)  can 
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be  rearranged  to  give  an  explicit  expression  for  the  algorithmic  yield: 

»•  =  (^)  •  (2) 

With  this  equation,  given  a.  b.  and  a  data  set  such  as  that  shown  in  Fig¬ 
ure  1,  obtained  from  CORRTEX.^'^l  SLIFER,  or  another  time  of  arrival 
measurement  technique,  one  can  reformulate  the  data  as  algorithmic  yield 
vs  time,  as  is  shown  plotted  in  Figure  2.  Again,  the  quantities  are  scaled 
appropriately  using  the  official  yield  for  the  event.  There  is  an  extremely 
steep  region  of  the  curve  before  about  0.23  scaled  ms.  followed  by  a  rela¬ 
tively  flat  region  between  0.23  and  0.45  scaled  ms,  which  is  in  turn  followed 
by  an  almost  linear  rise  in  out  to  about  1  scaled  ms.  The  very  steep 
portion  of  the  curve  before  0.23  scaled  ms  is  not  valid  data.  This  data  is 
probably  due  to  energy  flow  along  channels  near  the  gauge  cable,  causing 
a  shock  wave  which  arrives  considerably  earlier  than  the  expected  free-field 
shock  and  which  is  of  sufficient  magnitude  to  crush  the  cable.  This  feature 
is  common  in  records  from  cables  located  in  or  near  the  device  emplacement 
hole.  The  rising  portion  of  the  curve  above  0.45  scaled  ms  is  beyond  the 
range  of  applicability  of  the  algorithm  owing  to  the  increasing  dominance 
of  medium  strength  effects  at  the  lower  shock  stress  levels.  Only  the  “flat” 
region  of  the  data  from  0.23  to  0.45  scaled  ms  is  used  in  the  determination 
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Algorithmic  Yield  vs  Time 


Figure  2:  The  ratio  of  aigoriihmic  yield  to  official  yield  vs  scaled  time  for 
the  same  event. 

of  the  hydrodynamic  yield,  The  yield  is  computed  by  arithmetically 

averaging  over  the  appro-ximaiely  constant  values  obtained  from  the  (Hj.f) 
pairs  within  the  algorithmic  region.i®^  although  least  squares  fitting  over 
the  interval  has  been  tried  es  well^®-  with  negligible  change  in  the  results 
obtained. 


The  method  described  here  leaves  room  for  a  certain  amount  of  inter- 
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pretation  by  the  individual  analyst,  particularly  in  the  determination  of 
the  time  interval  over  which  to  apply  the  appropriate  averaging  operation. 
Historically,  this  step  in  the  yield  determination  process  has  been  done  by- 
hand,  so  that  slight  differences  in  computed  yields  could  be  calculated  for 
the  same  event  by  different  analysts  working  from  the  same  data  set.  We 
have  developed  an  automated  procedure  for  choosing  the  so-called  “algorith¬ 
mic  window”  for  a  particular  data  set  such  as  that  shown  in  Figure  1  and 
have  implemented  the  procedure  in  a  comp'jter  program.  The  algorithmic 
time  windows  determined  by  this  procedure  for  the  43  events  considered  in 
this  report  are  very  similar  to  those  found  by  Filers  et  and  the  recal¬ 
culated  yields  using  these  new  windows  are  typically  within  1%  of  the  Filers 
et  al  values.  This  new  procedure  does  not  therefore  supply  any  advantage 
over  the  existing  one  in  terms  of  giving  better  answers,  but  simply  codifies 
the  process  of  selecting  the  appropriate  time  window,  thereby  removing  a 
measure  of  human  judgement  and  variability  from  the  yield  determination 
process,  in  a  manner  which  is  consistent  with  current  practice. 
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2  Analysis 

Our  primary  task  for  this  project  was  to  estimate  the  uncertainty  associ¬ 
ated  with  employing  theliydrodynamic  yield  algorithm  to  determine  nuclear 
yields  for  underground  events.  For  our  hypothetical  test  we  are  to  assume 
the  data  is  obtained  from  CORRTEX  measurements  from  cables  emplaced 
in  a  satellite  hole  configuration.  With  properly  placed  satellite  holes,  con¬ 
tamination  of  the  data  by  such  phenomena  as  pipe  flow  should  not  be  a 
problem  For  this  analysis,  we  assume  spherical  symmetry  for  the  outgoing 
shock  front,  and  that  the  event  is  fully  coupled. 

In  order  to  ascertain  the  overall  uncertainty  in  the  application  of  the  al¬ 
gorithm,  we  begin  by  estimating  the  primary  sources  of  uncertainty  and.-their 
magnitudes.  The  algorithm  is  based  on  a  very  simple  equation.  Eqn.  (2). 
with  only  two  independent  variables,  i?,  and  t,  and  a  very  simple  hypoth¬ 
esis,  namely  that  Eqn.  (2)  holds  for  some  range  of  t,  given  particular  fi.xed 
values  of  a  and  6,  independent  of  geologic  material.  Thus,  our  first  task 
will  be  to  obtain  uncertainty  estimates  for  the  measured  quantities  Rj  and 
i,  and  to  judge  the  degree  of  generality  with  which  the  above  hypothesis 
fits  the  observed  data  for  the  range  of  different  media  represented  in  the 
records  available  to  us.  These  estimates  will  be  based  largely  on  a  statistical 
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analysis  of  the  (/Zj.t)  data  in  our  database,  as  well  as  current  information 
on  the  stale  of  the  art  in  drilling  and  surveying  (for  determining  relative  lo¬ 
cations  of  gauges  and  the  explosion  center)  and  in  timing  technology.  After 
these  constituent  uncertainties  have  been  estimated,  the  three-point  analysis 
method  of  Binninger  and  WrighJ^‘^^1  will  be  applied  to  obtain  the  resulting 
overall  uncertainty  in  yiejd  determinations  produced  by  the  strong-shock  al¬ 
gorithm  The  three-point  methodology  is  an  analysis  tool  which  has  proven 
useful  in  evaluating  the  effects  of  the  individual  uncertainties  involved  in 
the  application  of  a  more  complicated  procedure.  We  will  also  attempt  to 
determine  whether  the  '"blackout  region”,  which  severely  reduces  the  like¬ 
lihood  of  getting  useful  shock  time-of-arrival  measurements  in  the  extreme 
radiation  environment  near  the  device,  places  any  significant  restrictions  on 
the  utility  of  the  strong-shock  algorithm  over  any  substantial  portion  of  the 
yield  domain  of  interest. 

There  are  a  number  of  factors  which  may  reduce  the  effectiveness  of  yield 
determination  techniques  based  on  the  strong-shock  algorithm.  First,  the 
time-of-arrival  (TCA)  measurements  must  be  taken  in  the  hydrodynamic  re¬ 
gion  of  the  blast,  a  region  whose  radius  scales  with  yield  and  may  therefore 
become  prohibitively  small  for  low-yield  events.  There  is  also  a  minimum 
stand-off  distance  inside  which  measurements  are  either  impractical  or  use- 
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less.  High  neutron  and  gamma  radiation  levels  close  to  -‘he  device  may 
cause  vaporization,  melting  or  ionization  of  gauge  elements.  .4lso,  in  the 
very  near-held,  non-ideal  effects  due  to  finite  source  size  (Eqns.  (1)  and  (2) 
assume  a  point  source)  and  finite  cavity  size  become  apparent.  These  are 
factors  which  even  in  the  best  of  circumstances  can  limit  the  usefulness  of 
algorithms  based  on  Eqn.  (1).  Uncertainties  in  the  relative  locations  of  the 
gauges  with  respect  to  the  effective  center  of  the  e.xplosion  may  also  influ¬ 
ence  our  ability  to  make  accurate  yield  determinations  using  this  algorithm. 
These  uncertainties  may  be  related  to  limits  on  the  absolute  accuracy  with 
which  the  gauges  may  be  placed  or  determined  (drilling  and  surveying), 
limited  knowledge  of  the  location  of  device,  or  offset  of  the  effective  cen¬ 
ter  of  the  explosion  (ECE)  from  the  location  of  the  original  device  center. 
It  is  conceivable  that,  under  the  circumstances  of  a  verification  situation, 
accurately  determined  point-of-origin  and  explosion-time  iiducials  may  not 
be  readily  available.  The  effect  of  each  of  these  elements  would  evidently- 
become  more  pronounced  for  small  yields.  Our  effort  on  this  project  has 
been  aimed  toward  quantifying  the  effects  of  each  of  these  and  other  factors 
on  our  ability  to  accurately  make  yield  determinations  from  hydrodynamic 
shock  location  measurements 

We  have  attempted  to  quantify  what  we  felt  were  the  primary  sources 
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of  uncertainty  associated  with  the  application  of  the  strong-shock  algorithm 
to  yield  determination  and  treaty  verification.  One  of  the  most  important 
questions  which  must  be  asked  pertains  to  the  primary  assumption  of  the 
strong  shock  algorithm  itself.  How  well  does  Eqn.  (1)  fit  data  for  hydro- 
dynamic  shock  location  in  a  wide  variety  of  different  geologic  media  using 
fixed  values  of  a  and  6,  and  with  what  degree  of  confidence  may  it  be  ap¬ 
plied  in  even  the  most  ideal  of  situations?  This  question  has  been  addressed 
elsewherej®’^^  but  our  approach  is  from  a  somewhat  different  perspective. 

The  parameters  a  and  b  have  both  been  empirically  determined.1^1  There 
are  a  number  of  considerations  involved  in  the  determination  of  the  parame¬ 
ter  6.  It  is  clear  that  the  location  of  the  shock  front  and  its  velocity  satisfy  an 
equation  like  (1)  at  every  instant  for  some  pair  of  values  o  and  6.  It  is  likely 
that  the  appropriate  value  for  6  will  lie  somewhere  between  0.4  (the  ideal  gas 
limit)  and  1  (the  linear  wave  limit),  with  the  lower  value  applying  near  the 
source,  and  6  —  1  as  /?  — *  oc.  Choosing  a  particular  value  for  b  effectively 
selects  a  portion  of  the  shock  location  history  corresponding  to  that  value. 
It  was  observedf^^  in  numerical  experiments  using  realistic  material  models 
that  the  value  6  =  0.475  is  sufficiently  small  to  select  an  algorithmic  region 
which  is  relatively  insensitive  to  geologic  media  (i.e.,  is  within  the  hydrody¬ 
namic  regime),  and  is  large  enough  to  avoid  many  of  the  perturbations  due 
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to  emplacement  geometry  (e.g.,  pipe  flow  or  cavity  effects)  which  may  be 
evident  in  the  near-source  region  for  many  tests.  This  value  has  been  used 
to  successfully  determine  yields  for  a  large  number  of  underground  nuclear 
tests  executed  in  a  range  of  different  geologic  media.^®’^  For  this  analysis, 
we  assume  a  fixed  value  of  6  =  0.475. 

2.1  Nuclear  shot  data. 

For  this  analysis,  we  have  considered  shock  time  of  arrival  data  from  43 
U.  S.  DOE  underground  nuclear  tests  occurring  between  1962  and  1983  and 
published  in  Filers  ei  Goldwire  and  Geill^l  and  Bass.l^  These  repre¬ 
sent  all  of  the  close-in  (hydrodynamic)  shock  time-of-arrival  records  readily 
available  to  us  for  use  in  this  analysis.  All  of  the  events  in  our  database 
had  yields  greater  than  forty  kilotons.  Most  of  the  data  were  obtained  us¬ 
ing  the  CORRTEX  experimental  technique,  while  the  remainder  come  from 
SLIFER  records.  !n  all  cases,  the  published  “official  yield’’  was  used  as  the 
actual  yield  for  each  event.  For  most  of  the  events  considered  here,  the  of¬ 
ficial  yields  were  determined  primarily  from  radiochemical  yields,  increased 
by  calculated  or  estimated  underground  enhancements.  However,  for  some 
events  hydrodynamic  yields  were  considered  in  the  official  yield  determina¬ 
tion,  so  that  an  entirely  independent  yield  measure  was  not  available  for 
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analysis. 

We  digitized  the  (W,t)  data  points  from  Goldwire  and  Geil^®!  for  use  in 
this  analysis.  From  these  data,  using  Eqn.  (1)  we  were  able  to  reconstruct 
what  we  felt  was  a  reasonable  representation  of  the  original  (/?,<)  data  sets. 
The  range  of  the  data  for  each  shot  to  be  used  in  our  analysis  of  the  standard 
yield  algorithm  was  taken  from  the  Filers  et  at  paper, so  that  we  would  be 
duplicating  as  closely  as  possible  the  the  “standard”  implementation  of  the 
algorithm.  This  range  of  useful  data  defines,  for  our  purposes,  the  so-called 
algorithmic  region  for  each  shot. 

With  this  information,  our  first  analysis  effort  was  to  attempt  to  place 
some  uncertainty  estimates  on  the  value  of  the  parameter  o  from  Eqn.  (1). 
With  aset  of  shock  radii  vs.  time  data,  Eqn.  (1)  may  be  solved  by  linear  least 
squares  for  which  minimizes  the  sum  of  the  squares  of  the  deviations 
between  the  official  and  computed  yield  at  each  data  point  in  the  algorithmic 
region.  Let  the  sum  of  the  squares  of  the  deviations  be  given  by 


then  it  is  easy  to  show  that 


o/it  — 


6/(1 -fr) 


(4) 


where  the  above  summations  are  taken  over  the  data  pairs  within 
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the  algorithmic  region  for  a  particular  shot.  Using  this  technique,  we  are 
able  to  find  a  best  value  for  the  parameter  a  (in  a  least-square  sense)  for 
each  of  the  43  events,  and  an  associated  standard  deviation  measuring 

the  spread  in  the  distribution  about  Oy,,  of  the  a,’s  generated  by  the  data, 
where 


From  these  results,  namely,  and  cry,,  for  each  shot,  we  may  determine 
the  best  a  for  all  of  the  events  together  by  constructing  the  weighted  average 
of  the  individual  values  according  to: 


where  the  On  and  an  are  now  simply  the  Oy,-,  and  o-y,,  from  Eqns.  (4)  and  (5) 
for  each  event,  and  the  summations  are  now  taken  over  the  43  shots.  The 
weighted  standard  deviation  a  which  is  to  be  used  as  our  measure  of  the 
uncertainty  associated  with  the  overall  parameter  a,  is  then  computed  from: 


2  _  E  (on  -  ay  fan 


In  a  similar  manner,  we  can  also  construct  the  best  a  and  associated  standard 
deviation  for  any  subset  of  the  database,  such  as  for  the  shots  taken  in 
alluvium,  basalt,  or  some  other  material.  Some  of  these  results  are  shown 
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Table  1:  Parameter  a  statistics  for  different  media. 


Category 

a 

a 

No.  cases 

Alluvium 

6.176 

0.043 

3 

Basalt 

6.162 

0.063 

2 

Granite 

6.221 

- 

1 

Rhyolite 

6.115 

0.094 

8 

Tuff: 

Partially  Saturated 

6.107 

0.054 

3 

Saturated 

6.064 

0.075 

25 

Tuff/Rhyolite 

6.143 

- 

1 

All  shots 

6.085 

0.082 

43 

Basalt  and  Granite 

6.164 

0.062 

3 

in  Table  1.  The  value  of  a  =  6.085  which  we  obtained  as  the  best  value  for 
all  43  events  agrees  well  with  the  Eilers  ei  value  of  a  =  6.093.  The  a 
associated  with  each  5  provides  a  measure  of  how  well  these  events  fit  within 
the  framework  of  the  algorithm  with  fixed  values  of  o  and  b. 

Note  that  this  estimate  of  the  uncertainty  in  a  necessarily  includes  a 
contribution  due  to  the  uncertainty  in  the  determination  of  Hq.  the  official 
yield.  For  example,  it  is  conceivable  that  a  single  value  of  a  is  in  fact  correct 
for  all  of  these  events,  and  that  the  observed  scatter  is  due  entirely  to  error  in 
the  official  yield  determination.  Although  this  is  indeed  possible,  the  stated 
official  yields  are  our  best  available  estimates  of  the  actual  yields,  and  in  the 
absence  of  other  independent  measurements,  these  must  be  regarded  as  the 
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true  yields  for  the  events  in  this  study. 

In  Figure  3  the  a/,j’s  computed  using  Eqn.  (4)  are  shown  plotted  as 
a  function  of  the  density  of  the  event  medium.  The  best  linear  fit  in  the 
least-squares  sense  to  these  points  is  also  shown,  suggesting  a  possible  small 
systematic  density  dependence  for  the  parameter  a. 

a  VS  Density 


s  _  _  _ ^ _ 

*  1.50  1.75  2.00  *.25  2.50  2.75  3.00 

Specific  Density 


Figure  3:  Parameter  a  plotted  as  a  function  of  medium  density. 


A  second  source  of  uncertainty  associated  with  application  of  the  hydro- 


dynamic  yield  algorithm  stems  from  uncertainty  in  our  knowledge  of 
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the  precise  location  of  the  shock  front  with  respect  to  the  effective  center 
of  the  explosion.  In  this  simplified  analysis  in  which  we  assume  that  the 
shock  is  spherically  symmetric,  the  primary  causes  of  uncertainty  in  R,  are 
the  uncertainties  in  the  locations  of  particular  points  on  the  CORRTEX 
or  SLIFER  cable  with  respect  to  the  device — or  more  properly,  the  effec¬ 
tive  center  of  explosion  (ECE) — due  to  survey  error,  and  uncertainty  in  the 
shock  front  location  due  to  cable  dynamic  crush  length  uncertainty.  Using 
state  of  the  art  technology,  drill  holes  can  be  surveyed  to  accuracies  of  less 
than  0.01%  of  the  hole  depth,!^^^  or  1  cm  at  100  meters.  Dynamic  crush 
length  uncertainty,  or  the  uncertainty  associated  with  the  exact  location  of 
the  crush  front  on  the  cable,  is  claimed  to  be  about  ±15  cm  for  SLiFER 
and  ±3  cm  for  CORRTEX.^^J  Our  analysis  of  the  data,  based  on  the  scat¬ 
tering  of  the  data  points  about  an  ideal  smooth  curve,  suggests  the  overall 
uncertainty  in  /?,  for  all  43  events  (including  both  SLIFER  and  CORRTEX 
records)  to  be  about  ±18.5  cm  at  the  2cr  level.  These  data  may  reflect  un¬ 
certainties  in  Rf  due  to  the  above  factors  as  well  as  other  factors  not  listed 
here,  such  as  nonspherical  shock  propagation  or  cavity  effects. 

In  a  verification  situation,  limited  knowledge  of  emplacement  geometry 
may  make  it  difficult  to  predict  a  priori  the  location  of  ECE,  and  accurate 
explosion  time  fiducials  (to  establish  1  =  0)  may  not  be  available.  These 
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problems  have  been  addressed  in  a  limited  sense  for  1-D  (spherical)  geometry 
by  Goldwire  It  seems  clear  that,  in  the  real  world  of  ic’  ''iple  dimensions, 
these  types  of  systematic  uncertainties  in  R,  and/or  !  ■••j-  . '.Ifectively  re¬ 
moved  by  employing  two  or  three  CORKTEX  cable  rui.  •!  .  'lOuted  around 
the  device  emplacement  hole.  Clean  records  from  two  indr:*i>ndent  measure¬ 
ments  could  be  used  to  establish  one  unknown  (either  •  ~  t  —  b  if  the 

other  were  known),  while  three  independent  measurements  could  establish 
both. 

Vv'e  estimate  that  the  uncertainty  in  the  measurement  of  t  itself  is  neg¬ 
ligible  compared  to  those  in  a  and  Rt  since,  with  modern  timing  devices, 
measurements  can  be  made  to  almost  au-bitrary  precision. 

2.2  Blackout  region 

We  also  investigated  the  possible  significance  of  the  so-called  “blackout  re¬ 
gion”  inside  of  which  TOA  measurements  would  be  precluded  because  of  the 
severe  radiation  environment  at  shot  time.  Heat  generated  by  the  intense 
neutron  and  gamma  ray  fluxes  close  to  the  radiating  nuejear  device  may  lead 
to  vaporization,  melting  or  ionization  of  the  TOA  gauge  elements,  prevents 
ing  valid  time  of  arrival  measurements  in  this  region.  Electromagnetic  pulse 
(EMP)  and  induced  Compton  currents  in  cables  and  dielectrics  can  also 
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severely  degrade  measurements  in  the  near-source  region.  VVe  obtained  es¬ 
timates  of  the  radiation  fluence  levels  at  which  CORRTEX  hardware  could 
sui  vive  and  successfully  operate,  typical  n-7  fluences  from  a  modern  ra¬ 
diative  device,  and  an  effective  absorption  coefficient  for  7-rays  in  a  dense 
rock  medium  (e.g  rranite).^^^  Ignoring  any  systematic  dependence  of  the 
source  spectrum  on  yield,  we  assume  that  the  output  fluence  from  the  device 
is  proportional  to  the  yield.  This  leads  to  the  equation, 

oc  W.  (8) 

where  rg  is  the  radius  of  the  blackout  region,  and  the  constant  of  propor¬ 
tionality  may  be  determined  from  the  above  estimates.  Carrying  out  this 
procedure,  we  obtain  a  curve  for  the  estimated  radius  of  the  blackout  region 
in  a  dense  rock  material.  This  curve  is  plotted  in  Figure  4.  Also  shown  are 
two  curves  showing  the  approximate  extent  of  the  algorithmic  window  for 
the  strong  shock  yield  algorithm,  using  0.16  ms/kt'^^  and  0.60  ms/kt*/^. 
This  (estimated)  inner  limit  of  the  algorithmic  region  lies  clearly  outside  of 
the  estimated  blackout  region,  suggesting  that  radiation  effects  of  this  sort 
probably  do  not  significantly  restrict  the  region  of  validity  of  the  strong- 
shock  algorithm  in  media  of  this  type. 
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Measurement  Region 


Figure  4:  Blackout  r&gion  and  algorithmic  window  (estimated). 

'i.3  TJhree-poiri'c  analysis 

We  'lOW  have  identified  certain  error  sources  associated  with  the  strong-shock 
yield  determination  methodology,  and  have  established  reasonable  estimates 
of  their  constituent  uncertainties.  In  the  simplified  situation  considered  here, 
the  primary  sources  of  uncertainty  are  believed  to  reside  in  the  measurement 
of  Rf  with  respect  to  the  effective  center  of  explosion,  and  in  the  degree  of 
universality  as.sociated  with  Eqn.  (1)  using  constant  values  for  a  and  b.  The 
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next  step  is  to  assess  how  these  constituent  uncertainties  affect  the  overall 
uncertainty  for  the  yield  determination  method  as  a  whole.  To  accomplish 
this,  we  employ  the  three-point  analysis  of  Binninger  and  Wright.  This 
methodology  provides  a  quick  and  easy  way  to  evaluate  and  rank  the  effects 
of  the  above  uncer*  xinties  which  are  involved  in  the  application  of  the  yield 
algorithm.  The  results  provided  by  the  three-point  method  have  been  shown 
to  compare  favorably  to  those  produced  by  more  mathematically  rigorous 
statistical  methods  such  as  Monte  Carlo.t^J 

To  implement  the  three-point  method,  we  first  produce  nominal,  upper 
and  lower  bound  values  (at  the  95%  confidence  level)  for  each  of  the  variables 
which  we  wish  to  include  in  the  uncertainty  analysis.  This  is  essentially  the 
2a  level  for  normally  distributed  random  variables,  so  here  we  use 
0=  {5.921,6.085.6.249} 


(9) 


Rs  =  {i?o  -  18.5cm.  Ro,  Ro  +  18.5cm}  . 


The  value  for  Ro  is  computed  using  Ro  =  with  a  =  6.085,  h  = 

0.475  and  t  =  0.38  ms  x  where  W  is  the  yield  in  kilotons.  This  value 
for  t  is  chosen  to  place  it  exactly  in  the  center  of  the  nominal  algorithmic 
window  which  ranges  from  0.16  to  C.60  scaled  milliseconds.  Clearly,  Ro  is 
yield  dependent,  while  a  is  not  :Fcr  a  given  yield  event,  we  now  construct 
a  table  of  the  yields  which  would  be  predicted  from  a  single  measurement 
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using  each  combination  of  values  in  (9).  As  an  example,  for  a  1  kt  shot, 
these  values  would,  according  to  Eqn.  (2),  generate  a  table  of  W  values  as 
shown  in  Table  2. 


Table  2:  One  kiloton  example. 


a 

values 


Ro  values 


3.658 

3.843 

4.028 

5.921 

0.8S2 

1.169 

1.529 

med  =  1.169 

6.085 

0.754. 

1.000 

1.308 

med  =  r.OOO 

6.249 

0.648 

0.859 

1.124 

med  =  0.859 

med  =  0./54  med  =  1.000  med  =  1.308 


To  quantify  the  effect  of  the  uncertainty  in  o,  for  example,  on  the  total 
dispersion  in  W,  the  variate  o  is  fixed  at  its  lower  bound  values,  and  then 
values  of  W  are  determined  according  to  the  values  of  Ro-  These  values  for 
W  are  ordered  according  to  ascending  values  of  /?©  as  shown  in  the  first  row 
of  the  matrix  in  Table  2.  Similarly,  the  values  of  W  are  deterrrjined  over  the 
range  of  Rq  for  a  fixed  at  its  nominal  value  (refer  to  the  middle  row),  and 
then  for  a  fixed  at  its  upper  value  (third  row). 

The  effect  of  o  on  the  dispersion  in  W  is  then  determined  from  the  sum  of 
the  absolute  differences  between  values  in  each  column  (measured  from  the 
median  in  each  column).  Observe  that  for  this  example  the  sum  of  absolute 
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differences  taken  over  the  three  columns  of  the  matrix  is  0.949: 


X>a  = 


=  10.882  -  0.7541  +  |0.754  -  0.754|  +  |0.648  -  0.754| 

+  11.169  -  1.0001  +  11-000  -  1.0001  -f  10.859  -  1.0001 


+  11.529  -  1.30S1  +  11-308  -  1.3081  -f  11.124  -  1.3081 


=  0.949. 


A  similar  procedure  is  used  for  /?,,  this  time  taking  the  absolute  differ¬ 
ences  across  each  row.  Summing  these  differences  for  each  row  in  Table  2 


we  get 


•Pfi. 


=  1.677. 

The  relative  contribution  of  the  uncertainty  in  a  or  Rg  to  the  total  disper¬ 
sion  in  IV  is  found  by  dividing  the  sum  of  the  absolute  differences  found  for 
each  component  source  of  uncertainty  by  the  total  sum  of  those  calculated 
for  both  a  and  72, .  Thus 


a  contributes: 


0.949-1-  1.677 


=  36%. 


R,  contributes: 


0.949 -hi  .677 


=  64%. 
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This  procedure  may  be  repeated  for  any  desired  nominal  yield  to  generate 
the  curves  in  Figure  5,  which  show  the  relative  contributions  of  the  a  and 
Rj  uncertainties  as  a  function  of  yield  for  0  <  IV  <  150kt. 


Uncertainty  Contribution  vs  Yield 


Figure  5:  Relative  contributions  of  uncertainties  in  a  and  R,  as  a.  function 
of  yield. 


To  estimate  the  uncertainty  in  yield  which  is  produced  by  these  sources, 
we  use  the  three-point  method  results  from  Table  2  to  construct  the  ap¬ 
propriate  log-normal  distribution  corresponding  to  these  results.  The  log- 
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normal  parameters  estimated  from  the  three-point  results  include  the  me¬ 
dian  W\i2  and  /?,  where  /?  is  a  measure  of  the  dispersion  exhibited  by  a 
log-normal  variate.  The  log-normal  distribution  f{W)  is  defined  by 

fiW)  =  W-^  exp{-(ln  W  -  In  .  ^>0.  (14) 

The  procedure  followed  to  quantify  these  parameters  is  as  follows: 


=  median  value  determined  from  array  of  W  values, 

,2  (15) 


r  = 


E.-i  {ln(W.i/VVV2))' 
N  -  1 


For  our  1  kt  example,  we  get  the  log-normal  parameters  —  1-0  and  /?  = 
0.273.  With  the  distribution  thus  defined,  we  can  compute  the  probability 
that  W  lies  below  a  certain  value,  say  Wp  ,  from  the  equations 


1 

FiWp)  =  exp.{-(lnr  -  In  Wi/2?/2l3^}  dx 


=  —f 

P  J- 


ln(W>/W,/,) 


exp{-x^/2P^}  dx 


(16) 


i  <  1  -f  erf 


fHWp/Wif^y 

V  n/2/?  > 


where  erf(2)  is  the  standard  error  /auction. Solving  Eqns.  (16)  for  F  = 
95%,  we  can  compute  the  95%  confidence  level  yield  threshold.  ^95%,  for 
a  given  yield  determination.  In  our  example,  if  a  1  kt  yield  determina¬ 
tion  were  made  by  applying  the  strong-shock  yield  algorithm  to  a  given  set 
of  TOA  data,  and  assuming  that  the  foregoing  assumptions  are  true,  then 
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solving  Eqn.  (16)  with  VV'i/j  =  1  kt  and  F  ;  95%  we  get  kV9s%  =  1-57  kt. 
That  is,  we  can  be  95%  certain  that  the  true  yield  of  the  measured  event  is 
below  1.57  kt.  A  similar  procedure  may  be  used  to  compute  the  any  confi¬ 
dence  limit,  for  any  nominal  yield  1^1/2-  Figure  6  shows  the  quantity 
(kV'p  -  W)/W  (the  fractional  yield  uncertainty)  plotted  as  a  function  of  the 
nominal  yield,  W,  at  the  95,  85,  75  and  65  per  cent  confidence  levels.  Each 
curve  shows  the  threshold  yield  (at  a  particular  level  of  confidence)  under 
which  we  would  expect  a  given  event  to  fall,  if  the  algorithmic  yield  had 
been  determined  as  W .  Thus,  for  example,  if  the  algorithmic  yield  for  a 
particular  event  is,  say,  100  kt  then  our  analysis  indicates  that  we  can  be 
95%  confident  that  the  true  yield  does  not  exceed  this  value  by  more  than 
27%.  For  large  yields,  the  fractional  uncertainty  at  the  95%  confidence  level 
appears  to  asymptote  to  about  24.5%.  This  is  because  the  uncertainty  in 
Ri  becomes  relatively  less  significant  for  large  yields  (see  Figure  5),  and 
the  yield  uncertainty  becomes  more  dominated  by  the  contribution  from 
the  algorithmic  parameter,  a.  which  is  independent  of  the  magnitude  of  the 
yield.  As  the  yield  approaches  zero,  the  uncertainty  in  /?,  becon-''*  relatively 
more  important  and  the  fractional  uncertainty  in  the  determination  of  IV  in¬ 
creases  dramatically,  with  fractional  yield  uncertainties  (at  95%  confidence) 
greater  than  37%  at  5  kt,  46%  at  2  kt  and  57%  at  1  kt.  For  lower  levels  of 
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confidence,  the  fractional  uncertainty  curves  are  commensurately  lower,  but 


still  show  a  sharp  rise  in  the  low-yield  limit. 


FVaclional  Yield  Uncertainty  vs  Yield 


Figure  6:  Fractional  yield  uncertainty  as  a  function  of  yield. 
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3  Conclusions 

Our  finding  of  yield  uncertainties  for  the  strong-shock  algorithm  in  the  range 
of  245-28%  for  events  with  yields  greater  than  about  40  kilotons  is  some¬ 
what  higher  than  the  overall  uncertainty  of  appraximately  20%  found  by 
Eilers  et  and  Goldwire  and  Geili^l  This  difference  is  likely  attributable 
in  large  part  to  the  differences  in  approach  taken  in  the  studies,  and  our 
limiting  value  of  24^%  for  high  yield  is  fairly  close  to  the  earlier  results. 
These  earlier  studies  also  do  not  explicitly  make  our  assumptions  of  a  fully 
coupled  event  and  spherical  shock  propagation.  However,  the  ranges  of 
yields  considered  effectively  insured  total  coupling,  while  knowledge  of  the 
emplacement  geometry,  accurate  time  fiducials,  and  appropriate  choice  of 
the  algorithmic  time  window  assured  that  for  the  range  of  data  considered 
the  shock  front  was  essentially  spherical  with  known  ECE.  In  a  validation 
environfhent  some  of  this  information  may  not  be  available,  so  steps  must 
be  taken  to  assure  that  measurements  in  the  hydrodynamic  (algorithmic) 
regime  are  also  in  the  region  of  spherical  shock  propagation,  and  a  suffi¬ 
cient  number  of  measurements  taken  to  be  able  to  pinpoint  the  ECE  and 
initiation  time. 


We  extend  our  results  into  the  limit  of  low  yields.  Here  we  find  that  as  the 
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yield  goes  to  zero,  the  fractional  yield  uncertainty  rises  sharply,  with  yield 
uncertainties  of  about  33%  at  10  let,  37%  at  5  kt,  46%  at  2  kt  and  57%  at  1  kt. 
Any  deleterious  effects  on  yield  analysis  due  to  possible  cavity  decoupling, 
timing  uncertainties,  asymmetric  device  placement  or  other  nonspherical 
geometry  would  of  course  be  magnified  at  low  yields,  increasing  the  level  of 
uncertainty  in  the  resultant  algorithmic  yields  in  this  regime. 

Our  analysis  shows  that  the  so-called  “blackout  region”  in  which  mea¬ 
surements  are  precluded  due  to  the  extreme  radiation  environment  does  not 
appear  to  significantly  restrict  the  application  of  the  algorithm  for  events 
surrounded  by  a  dense  geologic  medium.  If  a  significant  amount  of  the  mate¬ 
rial  between  the  device  and  the  measurement  regime  is  air  (as  in  the  case  of 
a  large  emplacement  cavity,  for  example),  the  range  of  this  exclusion  region 
would  of  course  be  extended,  possibly  precluding  valid  TOA  measurements 
in  the  algorithmic  region. 

For  fixed  b  =  0.475,  the  algc"ithmic  parameter  a  does  not  show  any 
strong  systematic  media  dependence,  although  there  is  e  "light  indication  of 
variation  with  material  density.  Thiis  finding  is  in  general  agreement  with 
that  of  Goldwire  and  Geil.^^^  Most  of  the  data  available  to  us  were  from 
events  executed  in  fully  saturated  tuff  (25  put  of  43  shots),  with  few  shots 
in  each  of  six  other  types  of  geology.  We  believe  that  there  is  sufficient  data 
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to  establish  that  a  does  not  vary  much  with  geology,  but  useful  conclusions 
about  the  magnitude  of  the  variation  which  might  be  used  to  fine-tune  the 
algorithm  require  more  data  in  media  other  than  tuff. 

Finally,  we  must  mention  that  our  analysis  is  based  upon  a  special  data 
set,  in  particular,  one  in  which  all  of  the  records  have  some  valid  adgorithmic 
data.  A  more  representative  sampling  of  data  records  may  include  records 
without  any  valid  data  in  the  algorithmic  region  (due  primarily  to  aspherical 
shock  propagation  or  cavity  effects).  Such  cases  may  be  difficult  to  detect 
without  independent  knowledge  of  the  device  yields.  Using  satellite  holes 
for  CORRTEX  reduces  this  risk  by  avoiding  pipe-flow  problems,  but  the 
effect  may  be  similar  for  asmalhyield  event  taking  place  in  a  relatively  large 
cavity,  since  the  algorithmic  region  is  necessarily  close  to  the  cavity  wall 
in  these  cases.  This  effect  may  be  increased  by  asymmetric  placement  or 
detonation  of  the  device.  Care  must  be  taken  to  apply  our  results  in  such  a 
way  that  our  assumptions  of  total  coupling  and  spherical  shock  propagation 
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